Read in Script

source("Safegraph_Script.R")
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion
rm(list=ls()[! ls() %in% c("rural_ipums","ruca", "subset_glass", 
                           "poverty_ipums", "pct_chg_aug", 
                           "industry_ipums", "aug_combined", 
                           "buffer1_vec","buffer2_perim_vec", 
                           "buffer2_vec", "cf_bf1", "cf_bf2", "gfperimvec",
                           "not_zcta", "rural_ruca", 
                           '%notin%', 'pp1', 'pp2', 'aug_combined', 'com')])

Read in Spatial Plots

plot.new() 
rasterImage(pp1,0,0,1,1)

plot.new() 
rasterImage(pp2,0,0,1,1)

Initially a single 20-acre brush fire, it rapidly grew and merged with two smaller fires that expanded to 11,000 acres during the night of September 27 into September 28.

The Glass Fire was fully contained on October 20, 2020, after burning over 67,484 acres and destroying 1,555 structures, including 308 homes and 343 commercial buildings in Napa County, as well as 334 homes in Sonoma County

Data Tables (Descriptives)

Fire Category Total ZCTA
GF 8
BF1 17
BF2 10

Using shape files from the GlassFire, 10 km buffers were used to create the first buffer (BF1), 10km surrounding BF1 creates the second buffer (BF2).

Each perimeter has unique ZCTA that is exclusive to their perimeter.

DT::datatable(com %>%
                filter(june != 'June') %>%
                mutate(new_date = as.POSIXct(strptime(date_range_start, "%Y-%m-%d"))) %>%
                mutate(week = as_date(new_date)) %>%
                group_by(city, zcta, week, fire_cat, 
                         Poverty_Dummy_75, Agr_Dummy_75,
                         Rural_Dummy_R) %>%
                summarize(sum_visit = sum(raw_visit_counts), .groups = 'drop'), 
              caption = 'Descriptive Statistics of August Mobility', 
              filter = "top",
              extensions = c("SearchPanes", "Select", "Buttons"),
              options = list(dom = "Btip",buttons = list("searchPanes")))
DT::datatable(com %>%
                filter(june == 'June') %>%
                mutate(new_date = as.POSIXct(strptime(date_range_start, "%Y-%m-%d"))) %>%
                mutate(week = as_date(new_date)) %>%
                group_by(city, zcta, week, fire_cat, 
                         Poverty_Dummy_75, Agr_Dummy_75,
                         Rural_Dummy_R)  %>%
              summarize(sum_visit = sum(raw_visit_counts), .groups = 'drop'), 
              caption = 'Descriptive Statistics of June Mobility', 
              filter = "top",
              extensions = c("SearchPanes", "Select", "Buttons"),
              options = list(dom = "Btip",buttons = list("searchPanes")))

Visualizations

Percent Change Using Fire Category and June as Baseline

pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta,fire_cat) %>%
  summarize(med = median(pct_chg_mon_f), .groups = 'drop') %>%
  ggplot(data =., aes(x = week, y = med, group = zcta)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=fire_cat), alpha = 0.2) +
  geom_point(aes(colour = fire_cat), alpha = 0.8) +
  ylim(-50,200) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +# week of Oct 12
  theme(axis.text.x = element_text(angle = 0)) + 
  labs(x="Date", y="Median Percent Change (Divide by 100)") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))  +
    annotate("rect", xmin = as.Date("2020-09-21"), xmax = as.Date("2020-10-12"), 
           ymin = -50.00, ymax = 150.00 ,
             alpha = 0, color= "orange") +
  ggtitle("Percent Change Mobility by Fire Category")
## Warning: Removed 33 row(s) containing missing values (geom_path).
## Warning: Removed 33 rows containing missing values (geom_point).

When looking at mobility: + Week of June 6, 2020 to the week of June 22, 2020 is used as a baseline + Percent change is calculated by week to week in each consecutive week in June + The median of percent change is then grouped by fire category in June and is then officially the baseline for August mobility + Then the percent change of August is then used for each consecutive week in August (August 31 - October 12) + A calculation of percent change is taken by using the weekly percentage change for August in relation to the fire category median percent change.

Fire Category Median Percent Change
GF 0.02
BF1 0.03
BF2 0.04

Percent Change of Mobility by Rural/ Fire Category and June as Baseline

Rural Category Total ZCTA
Rural 11
Urban 24 
pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta,Rural_Dummy_R, fire_cat) %>%
  summarize(med = median(pct_chg_mon_f), .groups = 'drop') %>%
  ggplot(data =., aes(x = week, y = med, group = zcta)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Rural_Dummy_R), alpha = 0.2) +
  geom_point(aes(colour = Rural_Dummy_R), alpha = 0.8) +
  ylim(-50,200) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +# week of Oct 12
  theme(axis.text.x = element_text(angle = 0)) + 
  labs(x="Date", y="Median Percent Change (Divide by 100)") + 
  facet_wrap(vars(fire_cat), nrow = 3) +
  ggtitle("Percent Change Mobility by Fire Category and Rural")
## Warning: Removed 33 row(s) containing missing values (geom_path).
## Warning: Removed 33 rows containing missing values (geom_point).

Using RUCA (Rural/Urban Commuting Codes) to classify what is urban/rural in a aggregated status.

https://depts.washington.edu/uwruca/ruca-uses.php https://www.ers.usda.gov/data-products/rural-urban-commuting-area-codes/documentation/

pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta, fire_cat) %>%
  group_by(Rural_Dummy_R, week, fire_cat) %>%
  summarize(med = median(pct_chg_mon_f), .groups ='drop') %>%
  ggplot(data =., aes(x = week, y = med, group = Rural_Dummy_R)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Rural_Dummy_R)) +
  geom_point(size = 2) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12")))  +# week of Oct 12
  theme(axis.text.x = element_text(angle = 0)) + facet_wrap(vars(fire_cat), nrow = 3) +ylim(-50,50) +
  ggtitle("Percent Change Mobility by Fire Category and Rural (Collapsed)") +
  labs(x="Date", y="Median Percent Change (Divide by 100)") 
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Warning: Removed 6 rows containing missing values (geom_point).

pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta) %>%
  group_by(Rural_Dummy_R, week, fire_cat) %>%
  filter(fire_cat == 'GF') %>%
  summarize(med = median(pct_chg_mon_f), .groups ='drop') %>%
  ggplot(data =., aes(x = week, y = med, group = Rural_Dummy_R)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Rural_Dummy_R)) +
  geom_point(size = 2) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +# week of Oct 12
  theme(axis.text.x = element_text(angle = 0)) + 
  labs(x="Date", y="Median Percent Change (Divide by 100)") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))  +
    annotate("rect", xmin = as.Date("2020-09-21"), xmax = as.Date("2020-10-12"), 
           ymin = -20, ymax = 60 ,
             alpha = 0, color= "orange") +
  ggtitle("Median Raw Visits by Date (GF/Collapsed)")
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).

Percent Change of Mobility by Agriculture/ Fire Category and June as Baseline

Agr Category Total ZCTA
Agr<75 26
Agr>75 9 
pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta,Agr_Dummy_75, fire_cat) %>%
  summarize(med = median(pct_chg_mon_f), .groups = 'drop') %>%
  ggplot(data =., aes(x = week, y = med, group = zcta)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Agr_Dummy_75), alpha = 0.2) +
  geom_point(aes(colour = Agr_Dummy_75), alpha = 0.8) +
  ylim(-50,200) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +# week of Oct 12
  theme(axis.text.x = element_text(angle = 0)) + 
  labs(x="Date", y="Median Percent Change (Divide by 100)") + 
  facet_wrap(vars(fire_cat), nrow = 3) + 
  ggtitle("Percent Change Mobility by Fire Category and Agr Industry")
## Warning: Removed 33 row(s) containing missing values (geom_path).
## Warning: Removed 33 rows containing missing values (geom_point).

Sum of ZCTA that have (Counts of under 0.5 and 0.5 - 0.99) / Total Count. This is coded as a new variable (Percentage Agg). Then to code which ZCTA in subset are in the upper quartile of agriculture - we split percentages up by quartile. Those that are in top 25% (5.725%) are dummy coded as upper quartile of Agriculture in subset.

# Agriculture BY Fire Cat
pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta, fire_cat) %>%
  group_by(Agr_Dummy_75, week, fire_cat) %>%
  summarize(med = median(pct_chg_mon_f), .groups ='drop') %>%
  ggplot(data =., aes(x = week, y = med, group = Agr_Dummy_75)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Agr_Dummy_75)) +
  geom_point(size = 2) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12")))  +# week of Oct 12

  theme(axis.text.x = element_text(angle = 0)) + facet_wrap(vars(fire_cat), nrow = 3) +ylim(-50,50) +
  ggtitle("Percent Change Mobility by Fire Category and Agr (Collapsed)") + 
  labs(x="Date", y="Median Percent Change (Divide by 100)") 
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Warning: Removed 6 rows containing missing values (geom_point).

pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta) %>%
  group_by(Agr_Dummy_75, week, fire_cat) %>%
  filter(fire_cat == 'GF') %>%
  summarize(med = median(pct_chg_mon_f), .groups ='drop') %>%
  ggplot(data =., aes(x = week, y = med, group = Agr_Dummy_75)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Agr_Dummy_75)) +
  geom_point(size = 2) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +# week of Oct 12
  theme(axis.text.x = element_text(angle = 0)) + 
  labs(x="Date", y="Median Percent Change (Divide by 100)") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))  +
    annotate("rect", xmin = as.Date("2020-09-21"), xmax = as.Date("2020-10-12"), 
           ymin = -20, ymax = 60 ,
             alpha = 0, color= "orange")+
  ggtitle("Median Raw Visits by Date (GF/Collapsed)")
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).

Percent Change of Mobility by Poverty/ Fire Category and June as Baseline

Pov Category Total ZCTA
Pov<75 26
Pov>75 9 
pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta,Poverty_Dummy_75, fire_cat) %>%
  summarize(med = median(pct_chg_mon_f), .groups = 'drop') %>%
  ggplot(data =., aes(x = week, y = med, group = zcta)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Poverty_Dummy_75), alpha = 0.2) +
  geom_point(aes(colour = Poverty_Dummy_75), alpha = 0.8) +
  ylim(-50,200) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +# week of Oct 12
  theme(axis.text.x = element_text(angle = 0)) + 
  labs(x="Date", y="Median Percent Change (Divide by 100)") + 
  facet_wrap(vars(fire_cat), nrow = 3) + 
  ggtitle("Percent Change Mobility by Fire Category and Poverty Status")
## Warning: Removed 33 row(s) containing missing values (geom_path).
## Warning: Removed 33 rows containing missing values (geom_point).

Sum of ZCTA that have (Counts of under 0.5 and 0.5 - 0.99) / Total Count

Find Quartile - Those that are in top 25% (11.415%) are dummy coded as upper quartile in subset.

# Poverty BY Fire Cat
pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta, fire_cat) %>%
  group_by(Poverty_Dummy_75, week, fire_cat) %>%
  summarize(med = median(pct_chg_mon_f), .groups ='drop') %>%
  ggplot(data =., aes(x = week, y = med, group = Poverty_Dummy_75)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Poverty_Dummy_75)) +
  geom_point(size = 2) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12")))  +# week of Oct 12

  theme(axis.text.x = element_text(angle = 0)) + 
  facet_wrap(vars(fire_cat), nrow = 3) + ylim(-75,100) +
  labs(x="Date", y="Median Percent Change (Divide by 100)") + 
  ggtitle("Percent Change Mobility by Fire Category and Pov (Collapsed)") 
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Warning: Removed 6 rows containing missing values (geom_point).

pct_chg_aug %>%
  na.omit() %>%
  group_by(week,zcta) %>%
  group_by(Poverty_Dummy_75, week, fire_cat) %>%
  filter(fire_cat == 'GF') %>%
  summarize(med = median(pct_chg_mon_f), .groups ='drop') %>%
  ggplot(data =., aes(x = week, y = med, group = Poverty_Dummy_75)) +
  scale_color_manual(values=c('Blue','Green', 'Red')) +
  geom_line(aes(colour=Poverty_Dummy_75)) +
  geom_point(size = 2) +
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +# week of Oct 12
  theme(axis.text.x = element_text(angle = 0)) + 
  labs(x="Date", y="Median Percent Change (Divide by 100)") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))  +
    annotate("rect", xmin = as.Date("2020-09-21"), xmax = as.Date("2020-10-12"), 
           ymin = -25, ymax = 60 ,
             alpha = 0, color= "orange") +
  ggtitle("Median Raw Visits by Date (GF/Collapsed)")
## Warning: Removed 2 row(s) containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).

Looking at PM25 Levels

Grouped by Rural/Urban

pct_chg_aug %>%
  na.omit() %>%
  inner_join(., subset_glass, by = c('zcta' = 'zcta', 'week' = 'date')) %>%
  select(zcta, fire_cat, week, PM25, wf_pm25_imp, bins, pct_chg_mon_z, Rural_Dummy_R) %>%
  group_by(week,zcta) %>%
  arrange(zcta) %>%
  ggplot(data = ., aes(x = week, y = PM25, 
                               color = Rural_Dummy_R)) +
  geom_point(aes(color = Rural_Dummy_R), alpha = .5) +
  scale_color_manual(values=c('Blue','Green')) +
  geom_smooth(se = F) + 
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +
  labs(title="PM25 Values by Date and Rural Status") +
  labs(x="Date", y="PM25 Values") +
  labs(color='Nation')  +
  annotate("rect", xmin = as.Date("2020-09-21"), xmax = as.Date("2020-10-12"), 
           ymin = 0, ymax = 100 ,
             alpha = 0, color= "orange")  + facet_wrap(vars(fire_cat), nrow = 3)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 34 rows containing non-finite values (stat_smooth).
## Warning: Removed 34 rows containing missing values (geom_point).

Grouped by Agricultural

pct_chg_aug %>%
  na.omit() %>%
  inner_join(., subset_glass, by = c('zcta' = 'zcta', 'week' = 'date')) %>%
  select(zcta, fire_cat, week, PM25, wf_pm25_imp, bins, pct_chg_mon_z, Agr_Dummy_75) %>%
  group_by(week,zcta) %>%
  arrange(zcta) %>%
  ggplot(data = ., aes(x = week, y = PM25, 
                               color = Agr_Dummy_75)) +
  geom_point(aes(color = Agr_Dummy_75), alpha = .5) +
  scale_color_manual(values=c('Blue','Green')) +
  geom_smooth(se = F) + 
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +
  labs(title="PM25 Values by Date and Industry") +
  labs(x="Date", y="PM25 Values") +
  labs(color='Nation')  +
  annotate("rect", xmin = as.Date("2020-09-21"), xmax = as.Date("2020-10-12"), 
           ymin = 0, ymax = 100 ,
             alpha = 0, color= "orange")  + facet_wrap(vars(fire_cat), nrow = 3)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 34 rows containing non-finite values (stat_smooth).
## Warning: Removed 34 rows containing missing values (geom_point).

Grouped by Poverty

pct_chg_aug %>%
  na.omit() %>%
  inner_join(., subset_glass, by = c('zcta' = 'zcta', 'week' = 'date')) %>%
  select(zcta, fire_cat, week, PM25, wf_pm25_imp, bins, pct_chg_mon_z, Poverty_Dummy_75) %>%
  group_by(week,zcta) %>%
  arrange(zcta) %>%
  ggplot(data = ., aes(x = week, y = PM25, 
                               color = Poverty_Dummy_75)) +
  geom_point(aes(color = Poverty_Dummy_75), alpha = .5) +
  scale_color_manual(values=c('Blue','Green')) +
  geom_smooth(se = F) + 
  scale_x_date(breaks = "1 week", limits = as.Date(c("2020-08-31", "2020-10-12"))) +
  labs(title="PM25 Values by Date and Poverty Status") +
  labs(x="Date", y="PM25 Values") +
  labs(color='Nation')  +
  annotate("rect", xmin = as.Date("2020-09-21"), xmax = as.Date("2020-10-12"), 
           ymin = 0, ymax = 100 ,
             alpha = 0, color= "orange")  + facet_wrap(vars(fire_cat), nrow = 3)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 34 rows containing non-finite values (stat_smooth).
## Warning: Removed 34 rows containing missing values (geom_point).